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We derive an analytic formula for the Young’s modulus in single-layer black phosphorus using 
the valence force held model. By analyzing the directional dependence for the Young’s modulus, we 
explore the third principle direction with direction angle 0t P = 0.268-7T besides armchair and zigzag 
directions. The maximum Young’s modulus value is in the third principle direction. More specifi¬ 
cally, the Young’s modulus is 52.2 Nm” 1 , 85.4 Nm , and 111.4 Nm" 1 in the armchair direction, 
zigzag direction, and the third principle direction, respectively. This new principle direction is of 
significance for future discussions of other anisotropic properties in the single-layer black phosphorus. 
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Few-layer black phosphorus (BP) is another interesting 
quasi two-dimensional system that has recently been ex¬ 
plored as an alternative electronic material to graphene, 
boron nitride, and the transition metal dichalcogenides 
for transistor applications^'— . This initial excitement 
surrounding BP is because unlike graphene, BP has a 
direct bandgap that is layer-dependent. Furthermore, 
BP also exhibits a carrier mobility that is larger than 
MoS^- The van der Waals effect in bulk BP was dis¬ 
cussed by Appalakondaiah et.aA First-principles calcu¬ 
lations show that single-layer BP (SLBP) has a band gap 
around 0.8 eV, and the band gap decreases with increas¬ 
ing thickness—£ For SLBP, the band gap can be manip¬ 
ulated via mechanical strain in the direction normal to 
the BP plane, where a semiconductor-metal transition 
was observed—'— 

The single-layer BP has a characteristic puckered 
structure, which leads to the two anisotropic in-plane 
directions. As a result of this puckered configu¬ 
ration, anisotropy has been found in various prop¬ 
erties for the single-layer BP, such as the optical 
properties,— — the electrical conductance,— the me¬ 
chanical properties^Ar— and the Poisson’s ratio i 16 ’ 18 : 19 

The present work focuses on the Young’s modulus of 
the SLBP. In all existing works, the Young’s modulus in 
the armchair direction is much less than that in the zigzag 
direction. For instance, the Young’s modulus from the ab 
initio calculations for the armchair and zigzag-directions 
is 28.9 Nm" 1 and 101.6 Nm” 1 in Ref [ 3 , or 21.9 Nm' 1 
and 56.3 Nm' 1 in Refjl5l or 19.5 Nm 1 and 78.0 Nm 1 
in Reffl9l 

From the above, in most existing works, the investi¬ 
gation of anisotropic properties is usually performed by 
comparing these properties in two principle directions, 
i.e., armchair and zigzag directions. However, in this 
work, we will disclose an additional principle direction in 
the SLBP, which will be referred to the third principle 
(TP) direction. The Young’s modulus has the maximum 
value in the TP direction, while the armchair and zigzag 


directions have the minimum Young’s modulus. 

In this paper, using the valence force field model 
(VFFM), we derive an analytic formula for the direc¬ 
tional dependence of the Young’s modulus in SLBP. Be¬ 
sides armchair and zigzag directions, we reveal the TP 
direction, in which a maximum Young’s modulus value 
is reached. 

The atomic configuration of the SLBP is shown in 
Fig. [l] The structure parameters were measured in the 
experiment— Two in-plane lattice constants are a\ = 
r 3 7 = 4.376 A and <12 = ?~24 = 3.314 A. The out- 
of-plane lattice constant is <23 = 10.478 A. The origin 
of the Cartesian coordinate system is in the middle of 
ri2- The x-axis is in the horizontal direction and the 
y-axis is in the vertical direction. There are four in¬ 
equivalent atoms in the unit cell di x a 2 of the SLBP, 
which will be chosen as atoms 1, 2, 3, and 6 in this work. 
The coordinate of these atoms are rq = (— uai,0, —ua 3 ), 
r*2 = (uai, 0, ua 3 ), r 3 = (0.5ai — uai, 0.5a2, ua 3 ), and 
r*6 = (—0.5ai + ttai, 0.5a2, —ua 3 ). The two dimension¬ 
less parameters are u = 0.0806 and v = 0.1017. The 
bond lengths from the experiment are d\ = r 2 s = ri§ = 
2.2449 A and d 2 = r 3 2 = 2.2340 A, and the two angles 
are Q\ = d 3 28 = 0.5357T and 0 2 = 6*321 = 0.56 77T. 

Several empirical potentials have been developed to de¬ 
scribe the atomic interaction for the SLBP, including the 
VFFM potential^ and the Stillinger-Weber potential— 
Both potentials were fitted to the phonon dispersion of 
the SLBP. The Stillinger-Weber potential includes some 
nonlinear properties, so it can be applied in molecular 
dynamics simulations of the SLBP. The VFFM is a lin¬ 
ear model, so it is suitable for the investigation of linear 
properties in the SLBP, like the elastic bending modulus 
studied in this work. The VFFM is convenient for de¬ 
riving analytic expressions for elastic properties thanks 
to its simplicity. An analytic expression is of help for an 
explicit understanding of the elastic properties. Hence, 
we will apply the VFFM to derive an analytic formula 
for the Young’s modulus of the SLBP. 
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TABLE I: Parameters (in eVA 2 ) for the VFFM potential from Ref [2~i1. 


K r 

K' r 

K e 

K' e 

K rr > 

Kr' 

K r g 

Ke 

Ke 

9.9715 

9.4598 

1.0764 

0.9341 

1.1057 

1.1057 

0.7207 

0.7207 

0.7207 




FIG. 1: (Color online) SLBP structure. There are three prin¬ 
ciple directions, i.e., armchair (blue arrows), zigzag (red ar¬ 
rows), and the TP (green arrows) directions. Color is with 
respective to the atomic z-coordinate. 


There are nine terms in the VFFM potential, 


K = 

\k (Adi) 2 ; 

(1) 

V' = 

v r 

\K(Ad 2 ) 2 -, 

( 2 ) 

Vg = 

X -K e d\ (Adi ) 2 ; 

(3) 

Vg = 

\Kgdid 2 (Ae 2 ) 2 ; 

(4) 

V rr / — 
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Kg = 
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(7) 

II 

7 ^K r g\Jdid 2 (Adi) (A 6 * 2 ) ; 

( 8 ) 

II 

\Ke v / d^(Ad 2 )(Ad 2 ). 
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term describes the energy associating with the variation 
of intra-group angles like 0234- The Vg term describes the 
energy variation due to the variation of the inter-group 
angles like 6*123. The V rr ' term describes the potential en¬ 
ergy for the simultaneous variation of two different intra¬ 
group bonds like r 23 and r 2 4 - The V^ r , term gives the 
potential energy for the simultaneous variation of bonds 
like r 2 i and r 2 3. The V r g term is for the energy associ¬ 
ation with the simultaneous variation of an intra-group 
bond like r 32 and an intra-group angle like $ 234 - The V'. s 
term gives the potential energy for the simultaneous vari¬ 
ation of an inter-group angle like 6*123 and an intra-group 
bond like r 23 . The V” e term gives the potential energy 
for the simultaneous variation of an inter-group angle like 
6*123 and an inter-group bond like ri2- All parameters are 
shown in Tab. Q] The unit of these parameters has been 
converted from dyne/cm in the original work to eVA -2 . 

To compute the Young’s modulus, the uniaxial uni¬ 
form strain e is applied to the structure. The magni¬ 
tude of the strain is e and the direction of the strain is 
e = (cos <fi, sin </>, 0). The angle </> is counted from the x- 
axis. Based on the VFFM potential, the strain energy 
density is 

W x So = 2^ (Ar 23 ) 2 + 2^ (Ar 28 ) 2 + 2^ (Ar 21 ) 2 

+ 4-^- (A6*32 8 ) 2 + 4-^ (A6* 3 2 i ) 2 

+ (A6*82i) 2 + 4 ” (Ar 23 Ar 2S ) 

+ ^~ 2 ~ (Ar 23 Ar 2 i) + (Ar 28 Ar 2 i) 

+ 4— K- (Ar 2 3A6*234) + 4—^— (Ar4 3 A6* 2 34) 

+ 4—^ (Ar23A6*i23) + 4 —(Ar 28 A6*i 28 ) 

+ 4^ (Ari 2 A0i2 3 ) + 4^ (An 2 A0 128 ) (10) 

where So = ai x a 2 is the area of the unit cell ai x a 2 . 
The right-hand side gives the total VFFM energy for a 
unit cell. 

The Young’s modulus can be obtained through its def¬ 
inition, 


The VFFM describes the energy variation of the system 
due to a small change in the bond length (A&,;) and the 
angle (Adi) with i = 1,2, which are induced by strain 
in the present work. The V r term describes the bond 
stretching energy for intra-group bond lengths like r 2 3 . 
The V' r term is the energy corresponding to the bond 
stretching for inter-group bond lengths like ri 2 . The Vg 


E 


d 2 W 

de 2 

^[2 K r (a 2 + a 2 ) + 2 K' r a 2 2 

DO 

4 Kg/3 2 + 4 Kg (/3 2 + /3 2 ) + 4AT rr /aia3 
AK' rr ,a 2 (ai + a 3 ) + 4 K r g (ai + a 3 ) /3i 
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+ 4 K' rS (otifo + a 3 / 3 , 4 ) + 4 K” g a 2 (p 2 + /?4)]> (11) 


where the strain-induced variations for the bond length 
and the angle have been expressed as linear functions of 
strain; i.e., Adi = c^e and A Oi = /3ie, with i = 1,2, 3,4. 
We have introduced six geometrical coefficients, 


a-i 

Pi 


dr 23l _<9r 12| _<9r 28| / 10 n 

e—0; 0^2 |^o; «3 - q_ |e—0; (12) 


de 

d0 3 28 | 

de 1 


de 


de 


,.o;A = %-|..n;/3, = %i|..o<13) 


de 


de 


The x-axis is rotated to the strain direction e. The 
coordinate for the vector, r = ( x,y,z ), in this new coor¬ 
dinate system is 


1 X 4, ' 


v+ 

\ z * / 

- 


— sin ej> cos <j> 0 

0 0 1 



(14) 


where the under script <j> denotes the new coordinate 
system with x-axis along the strain direction. In the 
new coordinate system, the strain is applied along the 
x-direction, so the effect of the strain is as follows 


( x e \ 

(l + e 

0 


II 

0 

1 

° 

\ z e ) 

l 0 

0 

1/ 


' x <!> \ 

V ^ > 

V H ) 


(15) 


where the under script e denotes the coordinate under 
strain. The first derivative of the vector is ^ = x^e^. 
As a result, the first derivative of the bond length is 

dr f I _ _ 1 _=> _ dfc_ _ 
de l e —0 r € e de r 

For the bond length d\ = r 23 , we have £230 = 
(0.5 — 2 u) ai cos <j> + 0.5a 2 sin^, so we get 

_ dr 23 , _ 1 2 

a l — ~5 le—0 — a '230 

de r 2 3 

1 9 

= — [(0.5 — 2u) ai coscj) + 0.5a 2 sin^i] . (16) 

«i 


For the bond length d 2 = r 2 h we have £210 = 

— 2 uai cos </>, so we get 

a -2 = = — xl H = (2uai cos <(>) 2 . (17) 

de r 2 1 v d 2 

After the strain is applied, it is obvious that 7-23 ^ r 2 s, so 
we have another bond length in the deformed SLBP, i.e., 
d 3 = r 28 . We have 228 4, = (0.5 — 2 u) a± cos 5a 2 sin <j>, 
which leads to 

_ dr 2S , _ 1 2 

a 3 ~ — e—0 — - x 284> 

de r 2 8 

1 9 

= — [(0.5 — 2u) a\ cos0 — 0 . 5 o 2 sin(/)] . (18) 

«i 


We consider the strain induced variation for the angle 
0i = d 3 28 ■ According to the definition, cos $328 = n 23 ■ 



FIG. 2: (Color online) Direction-dependent Young’s modu¬ 
lus for SLBP. The total Young’s modulus (black solid line) 
is mainly contributed by the V r potential term (blue doted 
line). The other eight potential terms (gray lines) contribute 
less than 10% to the total Young’s modulus. The maximum 
Young’s modulus is in the TP direction at <j> t p = 0.2687T. 


n 281 we have 


Pi 


d9 1 1 d 

IT = — : —H - o - cos ^ 1 
de sin de 

1 / 2x2 3 4>x 2 8ct> 2ai 

sin 61 \ d\ di 



Analogous derivation gives 


(19) 


P2 


dd 2 dd 32 \ 

de de 

1 2a; 2 i02 ; 230 


sin 0 2 


d\d 2 




( 20 ) 


We find that $321 yf d 3 2 i in the deformed SLBP. As a 
result, we have 



do 821 

de 

^X 2 l(j)X 2 8 (j) 


sin 62 


d\d 2 




( 21 ) 


Inserting the above geometrical coefficients into 
Eq. (fill) , we obtain the Young’s modulus for the SLBP. 
Fig. [2] shows the directional dependence for the Young’s 
modulus. The contribution to the Young’s modulus from 
all of the nine VFFM potential terms are displayed by 
individual curves. The Young’s modulus in the arm¬ 
chair direction is is much less than the Young’s mod¬ 
ulus in the zigzag direction. Similar anisotropy in the 
Young’s modulus has also been reported in several pre¬ 
vious studies ^— -16 ! 18 though the obtained values show 
variability between the different studies. The difference 
is probably due to different computational methods and 
potentials that have been used in different studies. In 
present work, the Young’s modulus is 52.2 Nm -1 in the 
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FIG. 3: (Color online) Illustration for cj> t p ^ . The cross¬ 

over between <J> = (j>(9 1 ) (blue solid line) from Eq. (1251) and 
<t> = ^2 ( re< f dashed line) happens at 9i = 0.5-7T. This is the 
only point satisfying <f > tp = 


armchair direction and 85.4 Nm ” 1 in the zigzag direc¬ 
tion, which is within the range of the previously reported 
values 

Fig.[2]shows that the Young’s modulus has a minimum 
value at <f> = 0 (armchair direction). However, this curve 
clearly demonstrates that the Young’s modulus actually 
has a minimum point at <j) = 7t/ 2 (zigzag direction), al¬ 
though the zigzag direction has larger Young’s modu¬ 
lus value than the armchair direction. The maximum 
Young’s modulus (111.4 Nm -1 ) is in the direction with 
</>t p = 0 . 2687 T. We call this additional principle direction 
as the TP direction. 

The Young’s modulus is mainly contributed by the V r 
term. We thus examined the Young’s modulus computed 
from the V r term as 


E r 


oo 

[(0.5 — 2u) ai ] 4 


1 


1 + tan </> 


^1 + 6 tan 2 ^ tan 2 </> + tan 4 tan 4 cj^j . (22) 

As a result, the extreme points are determined by the 
condition, = 0, which leads to 


tan(/> = 

0 ; 

(23) 

tan cj) = 

oo; 

(24) 

tan 2 cj) = 

1 — 3 tan 2 ^ 

(25) 

tan 2 ^ (tan 2 ^ — 3) 


These equations determine three principle directions in 
the SLBP. Eqs. (1231) and EH) give two local minimum 
for the Young’s modulus in the direction with = 0 
(armchair) or <f> = 7 t /2 (zigzag). Eq. E5l) determines the 
TP direction with </> tp = ± 0 . 2687 t, in which the Young’s 
modulus has the maximum value. 

We note that the TP direction is almost coincident 
with the bond direction like r *23 in Fig. [H i.e., </> tp ~ 
However, there is no guarantee for this equality accord¬ 
ing to Eq. (| 25 |) . This can be more explicitly illustrated 
in Fig. [3] which shows the functions (f> = 4>(0i) (blue 
solid line) from Eq. E5l) and (f> = ^ (red dashed line). 
The crossover between these two curves yields Q\ = 0 . 57 T, 
which is the only value satisfying 4> tp = t/. We ac¬ 
tually have = 0 . 5357 T for SLBP, which is not equal to 
0 . 57 T, but these two values are very close to each other. 
As a result, the TP direction is very close to the bond 
direction like F 23 . 

In conclusion, we have derived an analytic expression 
for directional-dependent Young’s modulus of SLBP. The 
Young’s modulus has the minimum value in the armchair 
direction, but the maximum Young’s modulus is not in 
the zigzag direction. Instead, we find the TP direction 
in the SLBP with <§> = 0.268-7T, along which the Young’s 
modulus is maximum. 
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